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Abstract 

We investigate the isochronous bifurcations of the straight-line librating orbit in 
the Henon-Heiles and related potentials. With increasing scaled energy e, they 
form a cascade of pitchfork bifurcations that cumulate at the critical saddle- 
point energy e = 1. The stable and unstable orbits created at these bifurcations 
appear in two sequences whose self-similar properties possess an analytical scaling 
behavior. Different from the standard Feigenbaum scenario in area preserving 
two-dimensional maps, here the scaling constants a and (5 corresponding to the 
two spatial directions are identical and equal to the root of the scaling constant 
5 that describes the geometric progression of bifurcation energies e n in the limit 
n — > 00. The value of d is given analytically in terms of the potential parameters. 

1. Introduction 

The present study arose in the context of applying Martin Gutzwiller's semiclassical trace 
formula PjJ |2] , and some of its extensions, to various model Hamiltonians and interacting 
fermion systems in the mean-field approximation, with the aim of describing prominent 
quantum shell effects semiclassically in terms of the leading classical periodic orbits with 
shortest periods. This approach was promoted by Strutinsky et al. 3j, who estimated 
the nuclear ground-state deformations from the shortest periodic orbits in an ellipsoidal 
billiard. Applications which involved the author were beats in the level density of the 
Henon-Heiles and related potentials jH E] , conductance oscillations in mesoscopic semi- 
conductor structures Ej , and the onset of mass asymmetry in nuclear fission [§. . We 
refer to the literature just quoted and to a recent monograph jH] for a detailed discussion 
of the extensions of Gutzwiller's theory that are adequate for treating the degenerate 
orbit families occurring in systems with continuous symmetries. Uniform approximations 
that become necessary in connection with symmetry breaking and bifurcations will be re- 
ferred to in the next section. The accumulated experience from these investigations is an 
astonishing performance of the periodic orbit theory in reproducing quantum-mechanical 
gross-shell structure, using just a few short periodic orbits in the appropriate semiclassical 
trace formulae. 

In the present paper, we will stay on a purely classical level and investigate the onset 
of chaos in the Henon-Heiles and related potentials through bifurcation cascades, the 
accompanying self-similarity of periodic orbits, and their analytical scaling behavior. 



1 



2. The Henon-Heiles potential 

We investigate here the role of the straight-line librating orbit A in the Henon-Heiles (HH) 
Hamiltonian 

H = \ (x 2 + y 2 ) + ~ (x 2 + y 2 )+e (x 2 y - i y 3 ) . (1) 

Introducing the scaled variables u = ex and v = ey, the scaled total energy e, in units of 
the saddle-point energy E* = l/6e 2 , becomes 



e = E/E* = 6 



^(w 2 + v 2 ) + V(u,v) 



3(w 2 + v 2 ) + 3(m 2 + v 2 ) + 6vu 2 -2v 3 . (2) 



The Newton equations of motion in u, v are 

u = -u(l + 2v), 

v = —v + v 2 — u 2 . (3) 

These equations, and therefore the classical dynamics of the HH potential, depend only 
on the scaled energy e as a single parameter. For our numerical investigations below, we 
have solved Eqs. (jSJ) numerically and determined the periodic orbits by a Newton-Raphson 
iteration using their stability matrix 

In the left part of Fig. we show the equipotential lines in the (u, v) plane. The lines 
for e = 1 intersect at the saddle points and form an equilateral triangle. Along the three 
symmetry axes (dashed lines) the potential is a cubic parabola as shown, e.g., along u = 
in the right-hand part of Fig. ^ The figure also shows the three shortest periodic orbits. 




Figure 1: The Henon-Heiles potential. Left side: equipotential lines in the (u,v) plane. 
The dashed lines are the symmetry axes. The three shortest periodic orbits A, B, and 
C (evaluated at e = 1) are shown by the heavy solid lines. Right side: scaled potential 
along u = 0. (After 



Henon and Heiles ^Oj have already observed that the classical motion in this potential 
is quasi-regular up to energies e ~ 0.5 and then becomes increasingly chaotic; when one 
reaches the saddle energy (e = 1), more than 95% of the phase space is covered ergodically. 
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The periodic orbits in the HH potential have been classified and investigated in detail by 
Churchill et al. [I2J and more recently by Davies et al. [13J. Up to e ~ 0.97 there exist only 
three types of periodic orbits with periods of the order of T = 2n (i.e., the fundamental 
period of the harmonic-oscillator potential reached in the limit e —>■ 0): the librations A 
and B, and the rotation C. Corresponding to the symmetry of the HH potential, orbits A 
and B occur in three orientations connected by rotations in the (u,v) plane about 27r/3 
and 47r/3. Orbit C maps unto itself under these rotations but has two opposite time 
orientations, whereas time reversal maps the orbits A and B onto themselves. The overall 
discrete degeneracies are thus three for orbits A and B, and two for orbit C. The orbits A 
are stable up to e ~ 0.81 where they become unstable at a period-doubling bifurcation. 
At higher energies they oscillate between stability and instability, undergoing an infinite 
number of isochronous bifurcations that cumulate at the saddle-point energy e = 1 where 
the period Ta becomes infinity. This bifurcation cascade is the main object of our present 
study. For e > 1 the orbits A do not exist any more, but a new set of three degenerate 
unstable librations across the saddle points come into being O HB1 Ej ; in accordance 
with Ref. ^1] we call them r. The orbits B are unstable at all energies; the orbits C are 
stable at low energies and undergo a period-doubling bifurcation at the energy e ~ 0.89 
beyond which they remain unstable. A, B, and C are the only generic periodic orbits 
in the HH system for e < 1; all other orbits are created through their bifurcations (and 
further sequential bifurcations). 




0.45 0.5 0.55 e 0.6 0.65 0.7 



Figure 2: Oscillating part of level density of the HH potential, Gaussian averaged over 
an energy range Ae = 0.0024, versus scaled energy e. Solid line: quantum-mechanical 
result (evaluated for e = 0.04). Dashed line: semiclassical result, obtained with the 
uniform trace formula of Ref. [2] including first and second repetitions of the three 
orbits A, B, and C. 
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In Ref. jl] the coarse-grained quantum level density of the HH potential was shown 
to exhibit a pronounced beating structure which could be reproduced quantitatively by 
Gutzwiller's semiclassical trace formula [1J including the orbits A, B, and C. In the limit 
e — > 0, the trace formula diverges due to the approaching harmonic-oscillator limit which is 
integrable and has SU(2) symmetry (and a continuous two-fold degeneracy of its periodic 
orbits). This divergence can be removed in a uniform approximation that has recently 
been derived for some specific cases of SU(2) and SO(3) symmetry breaking pj. In Fig. El 
we show a comparison of the oscillating part Sg(e) of the quantum- mechanical level density 
5g(e) of the HH potential with its semiclassical approximation obtained from the uniform 
trace formula of Ref. [5] including only the first and second repetitions of the three orbits 
A, B, and C. Both results have been coarse-grained by convoluting them with a Gaussian 
over an energy range Ae = 0.0024. The agreement is excellent up to e ~ 0.67. The errors 
at higher energies are expected to come mainly from the orbit bifurcations that have not 
been taken into account (and partly from inaccuracies in the quantum result [HE]). 

An attempt to include the bifurcations in the semiclassical approach has lead to our 
present study. Uniform approximations for isolated bifurcations of all generic types have 
been developed by Sieber and Schomerus [T^j and successfully applied to the semiclassical 
description of various systems with mixed dynamics (see also Ref. ^Hj for an alternative 
treatment of the three simplest bifurcation types). Interferences of two close- lying bifur- 
cations (so-called bifurcations of codimension two) were discussed in Ref. [T7j. However, 
none of these uniform approaches can be used in the present case of the orbit A, where 
an infinite number of bifurcations coalesce at the saddle energy e = 1. The uniform treat- 
ment of this bifurcation cascade in a semiclassical trace formula is a challenging task. We 
should mention that an infinity of orbit bifurcations in the integrable two-dimensional el- 
liptic billiard has recently been incorporated successfully into an analytical trace formula 
[TH] . The bifurcating orbit in the ellipse - the straight-line libration along the shorter di- 
ameter - and the orbit families created at its bifurcations are, however, of a rather simple 
nature, and it is not clear yet if we can apply the technique of Ref. [IE] to the present 
system. While working along this line [T^], it seemed worth while to investigate on a 
purely classical level the bifurcation sequence of the saddle orbit A in the HH and similar 
potentials, which bears a lot of resemblance to the famous Feigenbaum scenario observed 
in one-dimensional [23] and two-dimensional maps (2TJ [22] . We shall presently exhibit the 
self-similarity amongst the periodic orbits created at the successive bifurcations and show 
that it is quantitatively described by the same scaling constant 5 that accounts for the 
geometric progression of the bifurcation energies. Different from the Feigenbaum scenario, 
the constant S is given analytically here, but it is not universal in that it depends on the 
parameters of the potential. 

A convenient method to keep track of orbit bifurcations in a two-dimensional system 
is to plot the trace of their stability matrix, trM. Bifurcations occur whenever trM = 2. 
In Fig. El we show tr M for the primitive orbit A n (with its Maslov index n increasing by 
one unit at each bifurcation) and for the orbits born at its bifurcations. In the lowest 
panel, we see the uppermost 3% of the energy scale available for the orbit A. The first bi- 
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Figure 3: Trace of stability matrix M of orbit A and the orbits born at successive 
pitchfork bifurcations in the Henon-Heiles potential, versus scaled energy e. Subscripts 
are Maslov indices. From bottom to top: successively zoomed energy scale near e = 1. 



furcation occurs at e§ = 0.969309, where A5 becomes unstable (with trM^ > 2) and a 
new stable orbit R5 is born. At e$ = 0.986709, orbit A 6 becomes stable again and a new 
unstable orbit L 6 is born. In the middle panel, we have zoomed the uppermost 3% of the 
previous energy scale. Here the behavior of A repeats itself, with the new orbits R7 and 
Lg born at the next two bifurcations. Zooming with the same factor to the top panel, we 
see the birth of R g and L 10 . This can be repeated ad infinitum: each new figure will be a 
replica of the previous one, with all the Maslov indices increased by two units and with 
tr Ma oscillating forever. 

Note that we have only shown here the primitives (i.e., the first repetitions) of each 
orbit. The higher repetitions of A will also undergo regular bifurcations and exhibit a 
corresponding fractal behavior. (For instance, whenever tr M of an orbit becomes equal 
to —2, its second repetition will have trM = 2 and bifurcate.) This infinite proliferation 
of stable and unstable orbits creates an increasingly mixed phase space and paves the way 
to chaos, similarly to the well-known Feigenbaum scenario. 
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We should emphasize an important difference, though, to the Feigenbaum scenario of 
Refs. |2Dl ED E2j : the bifurcations investigated there were all period doublings. Following 
the new stable orbit born at each bifurcation to its next period-doubling bifurcation leads 
to the famous Feigenbaum tree with its fractal structure. In our present system, however, 
successive pitchfork bifurcations occur from one and the same orbit A. Due to the discrete 
symmetries of the HH potential, these bifurcations are not generic (which would imply 
period doubling) but they are isochronous (i.e., each new-born orbit at the bifurcation 
point has the same period as the parent orbit) j2H|. (We have tried to follow a sequence 
of period-doubling bifurcations in the HH system. However, this soon leads to very long 
orbits that become unstable very fast with increasing energy; their investigation after 
more than three doublings has turned out to be numerically very difficult. Also, the 
successive period-doubling bifurcations are not all of pitchfork type but seem to alternate 
between pitchfork and touch-and-go type.) 

One important result of the Feigenbaum theory was to establish the geometric progres- 
sion of the bifurcation values of the system parameter through a constant S that turned 
out to be universal for a certain class of maps. In the original work of Feigenbaum [20], a 
dissipative one-dimensional map was investigated. Successive studies in area preserving 
two-dimensional maps |21| 122] yielded a different value of 5. Translating to the present 
situation, we take the energy e as the system parameter and study the progression of 
the energy intervals 1 — e n between the cumulation point e = 1 and the n-th bifurcation. 
Note, however, that the new orbits are born in two different sets (see Fig. 01 below): the 
stable rotations R2 m -i with two time orientations, and the unstable librations L 2m that 
come in degenerate pairs lying symmetrically to the the saddle line containing the parent 
orbit (m > 3). It is thus necessary to study the corresponding sequences of bifurcations 
separately, so that we have to determine the ratios 

5n = (4) 

1 - e n+2 

separately for odd and even n, and to see if the values of 5 n become constant for large 
n. Averaging our numerical values 5 n from the range 7 < n < 12, we obtain 5 Q = 37.623 
from the odd n and 5 e = 37.633 from the even n. The standard deviation from their 
mean value 37.628 is 0.082, so that the difference between S e and 6 a is insignificant and 
we may conclude that the mean value 5 = 37.628 is unique. Below, we shall determine 
analytically an asymptotic value for S which confirms this numerical result. 

But let us first examine another important property of the Feigenbaum scenario for 
area preserving two-dimensional maps: two independent numerical constants a and /3 were 
found to describe the scaling of the fixed points at the bifurcations in the two directions 
of the map [22] • I n our Hamiltonian system, Poincare surfaces of section - chosen through 
the u or the v axis - would provide us with two-dimensional area preserving maps and 
allow us to study the evolution of the fixed points with energy. However, such a study is 
again hampered by numerical inaccuracies. Instead, we found a geometrical self-similarity 
of the periodic orbits born at the bifurcations that reflects the fractal pattern of the fixed 
points in the Poincare maps and can be analyzed numerically with higher accuracy. 
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Figure 4: Orbits bifurcated from the vertical A orbit in the Henon-Heiles potential at energy 
e = 1, shown with increasing Maslov indices from left to right. Dashed lines: rotations R,2 TO -i, 
solid lines: librations L/2 m (m > 3). Only one libration orbit is shown for each L2 m ; its partner 
is obtained by reflection at the vertical symmetry line u = 0. Top panels: successive scaling 
of u axis from left to right with the factor 0.163. Bottom panels: successive scaling of both 
axes with the same factor; along the v axis only the top part starting from v = 1 is shown. 



In Fig. 0] we show the shapes of the orbits born at the isochronous bifurcations of 
orbit A, with increasing Maslov indices from left to right; all were evaluated at the barrier 
energy e = 1. We chose them here to be oriented along the v axis on which their parent 
orbit A is lying. The closer they are born to e = 1, the smaller has their amplitude in the 
transverse u direction developed when they reach the barrier energy. Therefore, in the 
upper part of the figure, the u axis has been zoomed by a factor 0.163 from each panel 
to the next, in order to bring the shapes to the same scale. The orbits look practically 
identical in the lower 97% of their vertical range, but near the barrier (v = 1) they make 
one more oscillation in the u direction in each generation. In the lower part of the figure, 
we have zoomed also the v axis by the same factor from one panel to the next and plotted 
the top part of each orbit, starting from v — 1. In these blown- up scales, the tips of the 
orbits exhibit a perfect self-similarity. The fact that the same scaling factor was used in 
both directions means that we find the two scaling constants a and (3 to be identical here. 
We shall derive them analytically below and show how they are related to 5. 

Note that although the parent orbit A becomes non-compact and non-periodic for 
e > 1, all periodic orbits bifurcated from it survive up to arbitrary energy, becoming 
more and more unstable; in spite of the non-compactness of the Hamiltonian for 
e > 1 they stay in a finite region of space. At e = 1, all the orbits R.2m-i have become 
inverse-hyperbolically unstable, whereas the L 2m remain direct-hyperbolically unstable. 
Vieira and Ozorio de Almeida ^1] have also determined some of these orbits at e > 1, both 
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numerically and semi-analytically using Moser's converging normal forms near a harmonic 
saddle. They showed, like Davies et al. PHI , that with increasing energy e > 1 these orbits 
come closer and closer to the librating orbits r that oscillate across the saddles. (In Ref. 
[13], our orbits R5, L 6 , R 7 , and L 8 were named i a , %, i c , and id, respectively, and r was 
called S.) In Fig.EJwe present the orbits Rn and r at e = 1.00006625, projected onto four 
different planes in the phase space. The mixed projections (u,p u ) and (v,p v ) correspond 
exactly to the results shown in Ref. [T3]. Note in the upper left panel of the figure, how 
the Rn orbit (solid line) winds around the torus of the r orbit (dotted line). 
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Figure 5: Periodic orbits evaluated at the energy e = 1.00006625, projected onto four 
different planes in the phase space. Solid lines: orbit Rn, dotted lines: orbit r. Rightmost 
panels: zooming the v and p v axes near the saddle point. 



We shall now proceed to derive the analytical values of the scaling constants 5, a, and 
(3. We give here only the main idea of the derivation using intuitive arguments; a more 
detailed mathematical analysis will be presented in a forthcoming publication [21]. The 
key to the understanding of the geometric progression of the bifurcation energies e n is a 
plot of tr M not versus energy e but versus the period T. This is shown in Fig. |3 We see 
that for large T, the quantity trM of the orbit A exhibits a perfectly periodic sinusoidal 
dependence on Ta (which was noticed already in Ref. |E]). The asymptotic period of 
these oscillations is found here numerically to be AT = 3.6276 ± 0.0003. 

This result can be analytically derived by linearizing the equations of motion (|Tj) 
around the periodic orbit A which is a solution, e.g., with u^t) = 0. The one-dimensional 
motion of this orbit in the v direction is given by the scaled potential (cf. Fig. 

6V Q (v) = 3v 2 -2v 3 . (5) 
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Figure 6: Trace of the stability matrix M of the orbits A (heavy line), B, C, and the 
orbits R2m-i) L2m ( m > 3) born at successive pitchfork bifurcations of orbit A in the 
Henon-Heiles potential, plotted versus their individual periods T. AT is the asymptotic 
period of the curve trM^ (Ta) for large Ta- 



Hence we find from energy conservation 

Hv) = V3 === s - ( 6 ) 

Jvi \ / e — 3s z + 2s 6 

This integral can be expressed in terms of an elliptic integral of the first kind, F, as 

dv 



t(v) = y/3/2 = = ^6/(^-^(7,^). (7) 

Hereby (i = 1, 2, 3) are the zeros of the equation 6 V (vi) = e; Vi and v 2 are the turning 
points of the A orbit (see Fig. [TJ. The arguments of the elliptic integral are 

7 = arcsin y/(v - v{)/ (v 2 - v 1) , k 2 = (v 2 -Vi)/(v 3 - v x ) . (8) 

The period of the A orbit thus becomes 

T A = 2^/3/2 =2 v /6/(v 3 ~v l )K(k 2 ), (9) 

J vi V(f - Vx)(V2 - V){V 3 - V) 

where K is the complete elliptic integral of the first kind which diverges when the modulus 
k 2 becomes unity; this happens at e = 1 where v 2 — v s — 1. The function t(v) in (j2J) can 
be inverted in terms of the Jacobi elliptic function sn(s, k 2 ) to yield the exact solution for 
the v motion of the A orbit: 

v A (t)=v 1 + (v 2 -v l )sn 2 (s,k 2 ) 1 (10) 
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where s is the scaled time variable 



S = ty/(v 3 -V 1 )/6. (11) 

Linearizing Eq. (JTJ) around the solution u — 0, t> = yields the equations of 

motion for small perturbations 5u(t), 5v(t) around the A orbit: 

Su(t) + [l + 2v A (t)]5u{t) = 0, (12) 
5v(t) + [1 - 2v A (t)] Sv(t) = . (13) 

These equations are of harmonic-oscillator type with time-periodic frequencies, usually 
named after Hill [2H] who investigated them in connection with the lunar theory (cf. Ref. 
|2j). With the particular solution (|TU|) for v A (t), Eq. (jTJ| is actually a special case of the 
Lame equation (see Ref. [2E] for an extensive discussion of Hill's and related equations). 
It is from this equation that the stability of the A orbit towards small perturbations Su 
can be derived. Magnus and Winkler j2E] have given an iterative scheme to solve Hill's 
equation using the Fourier expansion of the time-periodic coefficient, which has been used 
to calculate the stability of a linear periodic orbit in the diamagnetic Kepler problem 
|27| |2"5] . We shall present the application of this procedure to the present Hamiltonian 
elsewhere [21], and just anticipate here that the result for tr is a series 

tr M A (T A ) = 2 cos {u ± T A ) + ... , (14) 

where the dots indicate correction terms coming from the higher Fourier components of 
v A (t) in IjlOj) . All corrections have the same period 2tt/uj± as the leading term in ()14j) 
and change therefore only the amplitude and the phase of the oscillations in trM^. The 
frequency u± is given by the constant term of the Fourier expansion, which equals the 
time average of the coefficient in ()12|) over the period T A ; in the limit e — > 1 it goes to a 
constant: 

<f ± = ([l + 2v A (t)]) TA — > 3. (15) 
From this, we immediately obtain the asymptotic period 

AT = 2ix/uj 1 _ = 2n/V3 = 3.6275987 ... (16) 

which is found approximately from the numerical curve tr (T A ) shown in Fig. EJ 

This result can be intuitively obtained by the following reasoning. Note that the orbit 
A spends most of its time near the saddle point where v — 1; this is the more true the 
closer the energy comes to e = 1. Replacing v A (t) by its saddle-point value, the coefficient 
in (fT2"|) becomes uj\ = 1 + 2 v A = 3. We are, in fact, just speaking in this lowest-order 
approximation of a harmonic oscillation transverse to the A orbit 

8u(t) = uq cos (u±t + <p) = uq cos (v^3 t + 0) , (17) 

with a constant frequency uj±_ given by the curvature of the HH potential at the saddle: 

col = d 2 V(u,v)/du\ =0v=1 = 3. (18) 
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Figure 7: Orbit Lio evaluated at e = 1. Right panel: shape of the orbit in the (u,v) 
plane. Left panels: numerical results for v(t) (above) and u(t) (below) versus time t (in 
steps of 0.01), shown by solid lines. Dashed line above: VA{t) of parent orbit A taken at 
the bifurcation energy eio- Dashed line below: approximation (|17|) with uq = 0.001017 
and the phase eft fitted to the numerical result. 



In this limit, the stability matrix of the A orbit has exactly the value trM^Tyi) = 
2 cos (uj±Ta) (see, e.g., Ref. [5], Appendix C.2.1), in agreement with the leading term of 
the expansion (|14jl. That 8u{t) is harmonic with frequency u± = v3 over most of the time 
period is clearly seen in the results presented in Fig. [?l where we show the numerical 
solutions for u(t) and v(t) of the orbit Li evaluated at the saddle energy e = 1. The 
transverse u motion is, indeed, fitted extremely well over most of the period by Eq. ()17j) 
shown by a dashed line; the value of u will be determined below. 

The periods T n at which the bifurcations occur are thus given asymptotically by 
T^m-i ~ a + m AT and T 2m ~ b + m AT (see Fig. [HJ), with AT given by Eq. (fTrjj) and 
some constants a, b. The bifurcation energies e n are now easily found by the asymptotic 
expansion of the period T& © near e = 1, i.e., near k 2 = 1, where it diverges like 

T A ~ In (4/Vl - k 2 ) + ... ~ In (432/e) + . . . , (19) 

with e = 1 — e. (This is easily derived from a Taylor expansion of the turning points Vi 
and the quantity k 2 in powers of e.) Hence we find, for both even and odd n, 

1 - e n ~ 432 e~ T " oc e ~ nAT/2 oc e ~ nn/uJl - . (20) 

We thus obtain with (j3J) the analytical value of the asymptotic energy scaling constant S 

S = e 27T/uj± = e 27T/V * = 37.622367 . . . (21) 

which confirms the approximate numerical result given after Eq. (J3J). 
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The spatial scaling constants a and (3 describing the self- similarity of the new periodic 
orbits are derived along the same lines. In the limit vji(t) = 1, the equation for u(t) 
is harmonic even in the full equations of motion and decouples from that for v(t). 
Furthermore, the equation for v(t) is of second order in u so that to lowest order in small 
u oscillations, v (t) is not changed at all. Consequently, as long as the amplitude u of the 
transverse motion u(t) of the bifurcated orbits remains small, their v motion is "frozen" 
at the bifurcation point and given by the solution wa(^) in (fTUj). evaluated at e = e n . The 
transverse motion u(t) of the new orbit therefore carries off all the extra energy e — e n 
available above the bifurcation. That the v motion of these orbits is frozen above the 
bifurcation energies is seen in Fig. El by the fact that they appear as almost vertical lines 
(for large enough n), which means that their periods are practically constant. Also, in 
Fig. [7]we see that v(t) of the orbit L 10 obtained numerically at e = 1 (solid line) is, indeed, 
identical to VAif) evaluated at eio (dashed line, hardly distinguishable from the solid line). 
As a consequence of the frozen v motion, the energy of the u motion available for the new 
orbit at the saddle-point energy is just 1 — e n . In the harmonic approximation (|17p. this 
energy is equal to 

1 -e n = 3{u 2 + u 2 ± u 2 ) = 3uj 2 ± u 2 = 9u 2 . (22) 
Hence we find that the transverse amplitude of the new orbit at the saddle is given by 

u = • (23) 

Using the numerical result (1 — eio) = 9.305 x 1CT 6 , the amplitude for the orbit Li 
becomes from the above relation uq = 0.001017 which is, indeed, the value that fits the 
numerical result for u(t) in Fig. Now, from Eq. (|2*3*|) we find that the ratio of the 
amplitudes u of two successive generations of orbits for n — > oo tends to 

V(l - e n )/(l - e n+2 ) — > a = V6 = e^. (24) 

This confirms the numerical scaling constant 1/a — 0.1630335 . . . approximately used for 
the u scaling in Fig. HJ The scaling constant (3 for the v direction, finally, is obtained from 
expanding the potential (jHJ around the barrier at v = 1: 

6V (v) = l-3(v-l) 2 + ... (25) 

Since the v motion of the nth orbit is frozen, its tip near the barrier is just the turning 
point i> 2 evaluated at e n . Equating the frozen energy e n with 6 Vofa) to leading order in 
(j2EJ), we find the distance 1 — v% of the tip from the saddle point to be 

l-«2= yJ{\-e n )lZ. (26) 

Its ratio for two successive generations of orbits thus goes with n — > oo to the same limit 
as that found in (}2~l"j) for the u scaling. Hence (3 = a, as found numerically in Fig. 0] 

This concludes the derivation of our main result: the analytical value f!21[) of 6, and 
the relation a = (3 = y/5. There are further interesting observations, to be interpreted 
in future work (21]. For large n the actions of the bifurcated orbits at the saddle energy 
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e = 1 all tend to the action of the orbit A which is known analytically j4J: Sa = 6/5. Also 
at e = 1, the values for trM of all the orbits R.2m-i tend to the same value —4.183 (see 
their intersecting lines in Fig. EJ), whereas those of the L 2m intersect at the value +8.183. 

Finally, we add a remark about the r orbits that exist only for e > 1. They are 
oscillations transverse to the saddles (see Fig. in particular the upper right panel). In 
the limit e — > +1 they become one- dimensional harmonic oscillations with the frequency 
uj± = y/3. Their period is, in this limit, equal to T T = 2tv/uj±_ which is identical to AT 
in ([16)1 . Their transverse motion feels a negative curvature (corresponding to the passage 
over the saddle parallel to the orbit A) with the value 

-uf { = d 2 V(u, v)/dv 2 \ u=0 ,v=x = -1 , (27) 

as is seen directly from Eq. ()13j) with va = 1- They are therefore unstable, and their 
stability matrix has the trace (cf. Ref. 0, Sect. 5.6.3) 

trM r = 2 Cosh ( W ||T T ) = 2 Cosh (2tt/n/3) . (28) 

One of the eigenvalues of M r is thus A = e 2n ^ with the Lyapounov exponent x = 27r/\/3 
(or cr = x/T T = ^|| = 1). A is here identical with the scaling constant 5 (12 lj) for the 
bifurcation energies e n < 1 at which the orbits Ra m -i and L 2m are born. This exhibits 
once more the intimate connection ^2 C3] between the two types of periodic orbits near 
the threshold e = 1, already displayed in Fig. 03 

3. Results for other two-dimensional potentials 

We have investigated various two-dimensional potentials with one or more harmonic sad- 
dles and orbits that oscillate along straight lines towards the saddles. In all cases, we could 
find the same type of bifurcation cascades and the same self-similarity of the new-born 
orbits obeying the relation a = (3 = y/5. We shall give three examples below. 

3.1. A simple integrable case 

A separable potential with a minimum and one saddle is obtained if one omits the term 
proportional to x 2 y of the HH potential (JIJ to get 

V{x,y)= l -{x 2 +y 2 )- l -ey\ (29) 

and otherwise proceeds exactly in the same way as in Sect. 2. The straight-line A orbit 
along the v axis here sees the same potential as the A orbit in the HH potential and 
thus has the same solution VA(t) given in Eq. (|TU|) . Since the potential separates in u 
and v, the system is integrable. Nevertheless, the A orbit undergoes an infinite cascade 
of bifurcations cumulating at the energy e = 1. The u motion is strictly harmonic with 
period T = 2ir, corresponding to uj± = 1, and the trace of the stability matrix of A is 
exactly tr (Ta) = 2 cos (T^). The bifurcations thus occur at the periods T n = 2mc 
with n = 1, 2, 3, . . . This leads to a scaling constant 5 = e 2n for the progression of 
the bifurcation energies e n . The new periodic orbits born at the bifurcations are here 
degenerate families with trM = 2 whose tips scale in u and v with the constant v5 = e 77 . 
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3.2. The Barbanis potential 

Omitting the term proportional to y 3 in the HH potential yields a potential that has been 
studied in 1966 by Barbanis [2"9*] : 

V(x,y) = ^(x 2 + y 2 )-eyx 2 . (30) 

This potential has a minimum at x — y — and two saddles at the energy E* = l/8e 2 . 
There are two straight-line periodic orbits A oscillating through the minimum towards 
the saddles. Scaling the coordinates with a factor e and rotating the coordinate system 
such that one of the saddle lines becomes the horizontal u axis, the potential is 

V(u, v) = \ (u 2 + v 2 ) -= (2 u 3 - 3 uv 2 + V2v 3 ). (31) 

2 3\/3 ^ ' 




0.864 0.866 0.8658 0.866 

u u 



Figure 8: Periodic orbits in the potential of Barbanis [55], evaluated at the saddle-point 
energy e = 1. Top panel left: primitive orbits A (in two positions), B, C, and H in the (u, v) 
plane, shown by solid lines. The dashed lines give the equipotential curves at e = 1. Top 
panel right: the first generation of orbits bifurcated from the horizontal A orbit, shown on 
the same scale. Lower panels: extreme right tips of four generations of bifurcated orbits. 
The position of the saddle is at u = v3/2 = 0.866025 which corresponds to the right 
margins. The scaling factor between successive generations is 0.1084 in both directions. 
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Fig.|Hlshows the primitive periodic orbits in this potential. In the top left panel, we see 
the equipotential curves for e = E/E* = 1 by dashed lines: they form a parabola and a 
straight line intersecting at the two saddles. The solid lines indicate the shortest periodic 
orbits evaluated here at e = 1. Note that the two straight saddle lines that contain the 
A orbits are no symmetry axes (in contrast to the HH potential). The potential (|3T|) 
has only one symmetry line which halves the angle between the saddle lines and contains 
a straight-line librating orbit H 4 that is unstable at all energies. There is one curved 
librating orbit B 4 that oscillates transverse to the symmetry axis and is also unstable at 
all energies (like the orbits B in the HH potential). C is a rotating orbit that remains 
stable up to e = 2.415. The orbits A have the same behavior as in the HH potential and 
create an infinite set of new orbits at isochronous pitchfork bifurcations. 

Since the potential is not symmetric about the saddle lines, the new libration orbits 
born at the bifurcations come in pairs with equal Maslov indices 2n but shapes, actions 
and stabilities that differ at energies e > e2 n ; we call them here L2 n and \j\ n - The new 
rotations R,2 n -i again have two time orientations and thus a discrete degeneracy of two. 
The shapes of these new orbits resemble much those in the HH potential (apart from the 
asymmetry of the unstable pairs). Their self-similarity is shown in the lower four panels 
of Fig. El where their tips close to the saddle point u = V3/2 are displayed for each 
generation. The panels from each generation to the next are scaled by a factor 0.1084 
both in u and v direction. The bifurcation energies e n are again found to cumulate in a 
geometric progression, with a numerical scaling constant 5 = 85.1 ± 0.5. 

The values of tr M of all the primitive orbits are shown in Fig. |U] as functions of their 
period T. We find again a periodic behavior of tr (Ta), here with period AT = 4.443. 

10 i . • i ■ ■ — n 1 ^ 1 -r- ^ 1 — r ^ 1 1 




Figure 9: Trace of stability matrix M of orbits A (heavy line), B, C, H, and the new orbits 
born at successive pitchfork bifurcations of orbit A in the Barbanis potential, plotted 
versus their individual periods T. 
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The analytical values of these constants are found exactly in the same way as for the 
HH potential. In fact, the one- dimensional barrier seen by the A orbits in the Barbanis 
potential is, apart from a trivial scaling factor v3/2, identical to that in the HH potential. 
The frequency of the transverse oscillations across the saddles here is found to be cu± = y/2, 
so that AT = \f2-K = 4.4428829 . . . , and the scaling constant 5 of this potential becomes 

5 = = 85.019695 . . . (32) 

The scaling constants for the self-similar tips of the periodic orbits are therefore 

a = (3 = V5 = 1/0.10845266... (33) 

Both these values confirm the numerically determined constants. 

3.3. The quartic Henon-Heiles potential 

In Refs. jSl EH] a quartic Henon-Heiles potential was investigated that has a four-fold 
discrete rotational symmetry and four saddle points. It is given, in scaled coordinates, by 

V(u, v) = ~ (u 2 + v 2 ) - ~ (m 4 + v 4 ) + ~ u 2 v 2 . (34) 

We refer to the above literature for a discussion of the shortest periodic orbits (which we 
have renamed here to simplify the notation). We show the equipotential curves and the 
orbits A, B, and C in the top left panel of Fig. The orbits A have the same behavior 
as those in the standard HH potential, except that they oscillate between two opposite 
saddle points. They undergo again an infinite series of isochronous pitchfork bifurcations. 
The difference to the standard HH case is that here the nature of the new orbits born 
at the bifurcations alternates between rotations R n and librations L n / also amongst the 
stable orbits (odd n) and the unstable orbits (even n), as is shown for the orbits R 5 , L 6 , 
L7, and Rs labeled explicitly in Fig. E3 As a consequence, the self-similarity of their tips 
near the saddles becomes apparent only over two generations. Correspondingly, each of 
the lower four panels contains the tips of two successive generations of orbits, and the 
scaling from one panel to the next is done with the factor a 2 = (3 2 = 5. The frequency 
of the small oscillations across the saddles here is u± = 2, and the period of trM^ (T4) 
becomes asymptotically AT = 27r/ 'lj± = ir. The expansion of Ta near the saddle energy 
gives here, with e = 1 — e, 

T A ~ V2 ln (64/e) + . . . (35) 
For the bifurcation energies e n we have thus asymptotically (note the extra factor \/2) 

l-e n ~ 64e" T " / A (36) 

With T n ~ mi/uj^ —mr/2, we find for the energy scaling constant 

5 = e w/V2 = 9.2206125... (37) 

Its inverse is 1/5 = l /a 2 = 0.108452664 . . . , i.e., the scaling factor used in Fig. DH1 
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Figure 10: Periodic orbits in the quartic Henon-Heiles (HH4) potential, evaluated at 
e = 1. Upper left panel: equipotential curves (thin lines; dashed for e = 1) and the 
primitive orbits A, B, and C (heavy lines). Upper right panel: the first two generations of 
bifurcated horizontal orbits R5, Lq (above) and L7, Rs (below; orbit L7 is dashed), drawn 
on the same vertical and horizontal scales as the orbits in the upper left panel. Lower four 
panels: Tips of the bifurcated orbits R5 - R20 near the left saddle (u = —1), plotted on 
increasingly zoomed scales. The scaling factor from each pair of generations to the next 
is 0.10845 in both directions. 



The numerical iteration of the bifurcation energies e n was easier in this potential. The 
average value of the 5 n given by Eq. (j3J) in the range 7 < n < 16 gives a mean value of 
5 = 9.2203 with a standard deviation of 0.02. The asymptotic value of AT = it is reached 
for the interval T19 — T i7 with an accuracy of 7 digits. 
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4. Summary and conclusions 

We have investigated cascades of isochronous pitchfork bifurcations of straight-line lib- 
rational orbits in the Henon-Heiles and similar two-dimensional potentials possessing one 
or more harmonic saddles. The bifurcation energies cumulate at the saddle-point energy 
(scaled to be e = 1); their geometric progression yields a scaling constant 5 that can be 
calculated analytically in terms of the potential parameters. The periodic orbits born at 
the successive bifurcations exhibit a self-similarity corresponding to scaling constants a 
and (3 in the two spatial directions that turn out to be identical and equal to a — (3 — y/5. 

This result applies to a whole class of Hamiltonian systems with harmonic saddles 
and straight-line librational orbits oscillating towards the saddle points. The trace of 
the stability matrix of these orbits as a funcion of their period Ta, trM^T^), oscillates 
with an asymptotic periodicity AT = 2tc/u)±, where u\ is the transverse curvature at 
the saddles, see Eq. (flSjl. Combining this with the asymptotic energy dependence Ta ~ 
d\sx[c/ (1 — e)], one finds the asymptotic values e n of the bifurcation energies. From those, 
the analytic form of the scaling constant 5 in terms of the parameters u± and d is found 
to be 5 = e 27T ^ ±d . The Lyapounov exponent of the new unstable orbits r librating across 
the saddles for e > 1 is, in the limit e — > +1, given by x — < 2'^uj\\/lj± (or o = x/T T = uin), 
where — cuy is the negative parallel curvature at the saddles, see Eq. (J2"?|) . 

The well-known Feigenbaum scenario in two-dimensional area preserving maps [21} 122] 
differs from that studied here in three respects. First, the bifurcations discussed there are 
successive period-doublings which form a fractal tree. Second, their constants 5, a, and 
(3 are all different from each other and only known numerically. Third, these constants 
appear to be universal for a whole class of quadratic maps, whereas in the present case 
the constant S depends explicitly on the parameters of the potential. We should bear in 
mind that the Poincare maps corresponding to the Hamiltonian systems studied here are, 
of course, no simple quadratic maps. 

An isochronous pitchfork bifurcation cascade of a linear periodic orbit has been found 
also in the diamagnetic Kepler problem j2BHHI]- This orbit does not approach a saddle 
point but becomes infinitely long in the limit E — > — where its bifurcations cumulate. 
The self-similarity of the periodic orbits born at the bifurcations has not been investigated; 
the progression of the bifurcation energies is easily found to yield the constant 5 = 1 (cf. 
Refs. |3*2"1 The same value 5 = 1 is also found for the bifurcations of the short 

diameter orbit in the ellipse billiard, which are m-uplings with m > 2 cumulating in the 
limit m — > oo at zero excentricity (cf. Ref. jTHj and Ref. j^j, problem 5.3). 

It will be interesting to study also the modifications arising in connection with non- 
harmonic saddles, or with curved librating orbits approaching a saddle. 
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Note added in proof: 

The transverse motion u(t) [or v(t) in Sects. 3.2 and 3.3] of the bifurcated orbits R n , L n 
is given, close enough to their bifurcation energies so that its amplitude remains small, 
by the periodic Lame functions 1 Ec™ (at) and Es™(at). Hereby m is the number of zeros 
that can be uniquely related to the Maslov index n, and p and a are fixed numbers that 
depend on the potential. We find that the trigonometric expansions of Ec™(z) and Es™(z) 
given in Erdelyi et al. reproduce our numerical results to a high degree of accuracy (see 
Ref. [24 for details). 

1 see, e.g., Higher Transcendental Functions, Vol. Ill, A. Erdelyi et al., eds. (McGraw-Hill, 
New York, 1955), Chapter XV. 
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